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DIRECT  NUMERICAL  SIMULATION  OF  AXISYMMETRIC  JETS 


Introduction 

Experimental  investigations  in  the  last  decade  have  shown  that  large  spanwise  co¬ 
herent  structures  dominate  the  entrainment  and  mixing  processes  in  shear  layers  [l]. 
Recently,  it  has  become  possible  to  study  these  structures  by  direct  numerical  simula¬ 
tion  of  their  large  scale  features.  These  simulations  are  an  important  alternative  and 
supplementary  tool  in  the  basic  research  of  the  properties  of  fluid  flow  transitioning 
to  turbulence.  Since  a  simulation  can  calculate  values  for  all  the  primary  flow  field 
properties  as  a  function  of  time,  statistical  information  can  be  obtained  about  the  sys¬ 
tem  through  spatial  and  temporal  averages.  In  addition,  parameters  of  the  flow  can  be 
easily  varied.  While  the  experimental  conditions  may  not  be  fully  controllable  in  the 
laboratory,  the  simulation  conditions  are. 

Numerical  studies  of  coherent  structures  have  used  spectral  [2],  vortex  dynam¬ 
ics  [3j,  and  finite-difference  [4-6]  techniques.  Numerical  studies  of  the  evolution  of 
flows  similar  to  those  seen  in  the  laboratory  experiments  have  considered  both  two- 
dimensional  planar  and  axisymmetric  shear  layers.  Previous  finite-difference  calcula¬ 
tions  have  modeled  either  temporally-developing  mixing  layers  [4),  where  it  is  assumed 
that  the  relevant  vortex  dynamics  takes  place  in  relatively  compact  space  regions,  or 
spatially  developing  layers  [5-7],  which  represent  the  more  realistic  situations  occurring 
in  the  laboratory. 

In  previous  work  we  have  performed  finite-difference,  compressible,  spatially- 
developing  simulations  of  planar  shear  flows,  with  the  objective  of  investigating  asym¬ 
metries  in  mixing  [6|  and  the  basic  mechanisms  involved  in  the  reinitiation  of  the  insta¬ 
bilities  [7].  Here  we  describe  finite-difference  calculations  of  the  evolution  of  the  Kelvin- 
Helmholtz  instability  for  a  spatially-evolving  compressible  axisymmetric  jet  emerging 
into  a  quiescent  background.  The  instabilities  sustain  themselves  through  a  feedback 
mechanism  in  the  flow.  The  evolution  and  merging  of  the  downstream  structures  af¬ 
fect  the  inflowing  material  upstream,  thus  triggering  the  growth  and  shedding  of  new 
vortices  [7].  In  addition,  we  study  the  effect  of  small,  random  pressure  fluctuations  at 
the  nozzle  orifice  on  the  growth  of  the  mixing  layer.  These  fluctuations  model  inflow 
perturbations  in  experimental  flows  arising  from  turbulence  and  boundary  layers  in  the 
nozzle. 


Manuscript  approved  December  10,  1986. 


The  Numerical  Model 


The  numerical  model  used  to  perform  the  simulations  solves  the  two-dimensional 
time-dependent  conservation  equations  for  mass,  momentum  and  energy  for  an  ideal 


d(pV) 


=  -V  -pW  -VP, 


=  -V  .  cV  -  V  ■  PV, 
dt 


where  e  =  P/{'r  -  1)  +  is  the  internal  energy,  and  V,  P,  p,  and  7,  are 

the  velocity,  pressure,  mass  density,  and  the  ratio  of  specific  heats.  The  equations  are 
solved  using  the  Flux-Corrected  Transport  (FCT)  algorithm  [8]  and  timestep-splitting 
techniques.  FCT  is  an  explicit,  fourth-order,  finite-difference  algorithm,  which  ensures 
that  all  conserved  quantities  remain  monotonic  and  positive.  FCT  modifies  the  linear 
properties  of  a  high  order  algorithm  by  adding  diffusion  during  convective  transport  to 
prevent  dispersive  ripples  from  arising.  The  added  diffusion  is  subtracted  out  appropri¬ 
ately  where  not  needed  in  an  anti-diffusion  phase  of  the  integration  cycle  to  maintain 
high  order  accuracy.  With  this  approach,  no  artificial  viscosity  is  required  to  stabilize 
the  algorithm. 

The  model  uses  inflow  and  outflow  boundary  conditions  which  have  been  developed 
and  tested  for  these  types  of  multidimensional  FCT  calculations  [6,9] .  The  conditions 
ensure  the  proper  behavior  of  the  fluid  near  the  boundaries.  The  inflow  boundary 
conditions  specify  the  inflow  density  and  velocity,  and  use  a  zero  slope  condition  on  the 
pressure  at  the  inflow  boundary  to  define  the  energy  at  the  guard  cells: 


Pg  —  Pin  f  low  1 


^g  —  ^inflowt 


Pg  =  Pu 


where  Pi  is  the  pressure  at  the  first  (inflow)  cell.  This  adlows  pressure  differences  be¬ 
tween  the  jet  and  the  surrounding  fluid  to  generate  transverse  flows.  In  addition,  a 
short  inflow  plenum  is  modelled  by  including  a  portion  of  the  nozzle  within  the  com¬ 
putational  domain.  The  outflow  boundary  conditions  define  the  density  and  velocity 
at  the  guard  cells  by  means  of  extrapolations  from  the  last  two  cells: 


Pg  —  Pny 


Vg  =  2Vn  -  Wn-l- 


The  guard  cell  pressure,  in  turn,  is  defined  by  interpolating  between  the  boundary 
pressure  value  and  ambient  pressure  (assumed  at  infinite  distance  from  the  trsdling 
edge  of  the  nozzle): 

Pg  =  Pn+  -  Pn),  (5c) 
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where  V  is  either  the  radial  or  axial  coordinate,  and  the  subscript  j  refers  to  the 
trailing  edge  of  the  nozzle.  In  this  way,  the  outfiow  boundary  conditions  impose  a 
slow  relaxation  of  the  pressure  towards  the  known  ambient  value.  These  boundary 
conditions  allow  feedback  to  occur  between  the  downstream  vortices  and  the  inflowing 
material.  This  physically  realistic  approach  avoids  the  need  to  constantly  drive  the 
instability,  allowing  for  it  to  evolve  naturally  in  the  calculation.  No  subgrid  turbulence 
model  has  been  included.  The  simulations  are  expected  to  be  adequate  in  describing 
large  scale  features  of  gas  phase  flow  for  large  Reynolds  numbers. 

The  computational  grid  was  set  up  initially  and  held  fixed  in  time.  The  timesteps 
were  chosen  to  satisfy  the  Courant  condition.  The  basic  finite-difference  grid  used  120 
cells  in  the  cross-stream  (radial)  direction  and  220  in  the  streamwise  (axial)  direction, 
with  the  mesh  spacings  varying  in  the  ranges  0.05  <  AZ  <  0.52  cm  and  0.02  <  AR  < 
0.067  cm.  Figure  1  shows  a  schematic  diagram  of  the  grid.  The  cells  are  closely 
spaced  in  the  radial  (R)  direction  across  the  shear  layer,  where  the  large  structures 
form,  and  they  become  farther  apart  as  the  distance  from  the  shear  layer  increases 
for  R  >  Ro-  The  cell  separations  in  the  streamwise  direction  {Z)  also  increase  in  size 
as  we  move  downstream  and  upstream  from  the  trailing  edge  of  the  nozzle,  located 
a.t  R  =  Ro  =  0.9  cm,  Z  =  Zo  ^  1.67  cm.  This  takes  advantage  of  the  fact  that  the 
structures  merge  and  grow  downstream,  so  that  fewer  cells  are  necessary  to  keep  the 
same  relative  resolution  as  obtained  near  the  nozzle.  Convergence  of  the  results  was 
verified  by  checking  their  consistency  with  results  obtained  on  grids  with  440  x  120 
and  220  x  240  computational  cells,  which  doubled  the  resolution  in  the  axial  and  radial 
directions,  respectively,  relative  to  the  220  x  120  grid. 


rhe  Unforced  Axisvmmetric 


The  system  studied  is  a  high  speed  jet  containing  a  mixture  of  molecular  hydro¬ 
gen  and  nitrogen  emerging  into  a  quiescent  background  mixture  of  molecular  oxygen 
and  nitrogen,  with  a  jet  to  background  density  ratio  of  0.67.  The  system  was  initial¬ 
ized  with  a  step  function  axial  velocity  profile  at  a  uniform  temperature  (298  K)  and 
pressure  (1  atm).  The  jet  was  subsonic,  with  Mach  number  0.57,  and  corresponding 
velocity  2.0  x  10^  cm/s.  Figure  2  is  a  schematic  diagram  of  the  flow  configuration. 
The  perturbation,  which  initiated  the  instabilities  and  occurred  only  at  the  very  be¬ 
ginning  of  the  calculation,  is  a  very  small  cross-stream  pressure  gradient  generating 
vorticity  at  the  shear  layer  just  ahead  of  the  nozzle  edge.  This  disturbance  moves 
downstream  as  the  integration  proceeds,  generating  the  transverse  flows  which  trigger 
the  Kelvin- Helmholtz  instability.  Previously,  initial  sinusoidal  perturbations  along  the 
shear  interface  were  considered  [9].  The  current  approach  to  initiating  the  instability 
was  used  in  the  simulation  of  planar  shear  flows  [6,7],  and  is  closely  anaJogous  to  using 
a  delta  function  perturbation  at  the  center  in  an  idealized  periodic  simulation  involving 
two  equal  and  opposite  streams. 


Typical  features  of  the  flow  at  the  early  stages  are  shown  in  the  sequence  of  isovor- 
ticity  contours  in  Fig.  3.  Vortex  rings  develop  because  of  the  non-linear  growth  of  the 
instability.  This  occurs  at  an  essentially  fixed  distance  from  the  nozzle  edge,  somewhat 
less  than  one  diameter  D  =  2Ro,  at  z  =  (Z  -  Z^JD  a;  0.4.  The  newly  formed  struc¬ 
tures  move  along  the  interface,  interact  with  each  other  and  thereby  spread  the  vorticity 


until  the  center,  potential  core  region  disappears,  at  approximately  2  =  3.  The  struc¬ 
tures  are  displaced  vertically  by  a  low-frequency  modulation  of  the  shear  layer.  When 
appropriately  phased  with  the  structures,  the  modulation  tends  to  favor  the  merging 
of  three  structures  in  half  of  the  modulation  cycle,  where  they  will  coalesce  into  larger 
structures.  Since  the  jet  is  unforced,  the  low-frequency  modulation  can  only  be  due 
to  the  pressure  field  induced  by  the  larger  structures  downstream,  as  they  pass  near 
the  end  of  the  potential  core  of  the  jet.  As  the  effect  of  the  downstream  events  on  the 
inflowing  fluid  becomes  important,  there  is  an  interaction  between  the  basic  instability 
mechanism  at  the  shear  layer  and  the  feedback  mechanism.  The  flow  pattern  becomes 
subsequently  different  from  that  at  the  initial  stages  in  the  flow  development,  which  is 
dominated  by  the  shear  layer  instability.  A  striking  feature,  on  the  last  panel  of  Fig.  3, 
is  the  spatial  coherence  between  the  structures.  This  is  characteristic  of  the  feedback 
phenomena  between  the  downstream  events  and  the  inflowing  fluid  in  jets  [lO’j.  The 
distance  between  the  second  and  the  third  merging  locations  increases,  relative  to  the 
distance  between  the  second  and  the  first  mergings,  by  a  factor  of  the  order  of  2  -  3,  in 
good  agreement  with  the  results  of  jet  experiments  with  low-level  acoustical  excitation 
:i2i. 

Frames  at  much  later  times  in  the  development  of  the  flow  can  be  seen  in  Fig.  4.  A 
noticeable  feature  of  this  sequence  of  frames  is  a  regular  repeating  spatial  pattern  in  the 
evolution  of  the  flow.  We  have  performed  a  spectral  analysis  of  the  axial  (streamwise) 
velocity  fluctuations,  based  on  the  results  of  the  calculations  for  these  later  times. 
In  the  sequence  of  panels  in  Fig.  5,  we  show  the  results  of  the  analysis  at  various 
relevant  axial  locations  on  the  center  of  the  shear  layer  {R  =  Ro),  and  in  terms  of 
the  normalized  (Strouhal)  frequency  St  =  2ff/tfo./vo,  where  Oo  and  Va  are  the  initial 
momentum  thickness  of  the  shear  layer  and  the  jet  velocity,  respectively.  The  latter 
thickness  was  taken  to  be  effectively  defined  by  0o  =  Ai?  =  0.02  cm,  where  AR  is  the 
size  of  a  radial  cell  in  the  region  of  the  shear  layer  (within  which  the  step  of  the  initial 
velocity  profile  is  defined). 

Figure  5a,  at  z  =  0.4,  corresponding  to  a  location  near  the  nozzle  exit,  shows 
peaks  located  at  the  frequencies  Sto  as  0.100,  Sti  ss  0.050,  St2  as  0.025,  Sf+  as  0.125, 
and  St-  ss  0.075.  Here,  Stg  can  be  associated  with  the  shear  layer  instability  fre¬ 
quency,  in  good  quantitative  agreement  with  the  predictions  of  linear  inviscid  in¬ 
stability  theory  [13].  In  addition,  the  modulation  of  the  initial  shear  layer  with 
Sti  ss  Stg/2  and  Stj  ^  Stg/A  is  associated  with  the  first  and  second  subharmonic, 
while  St-i^ss  result  from  the  non-linear  interaction  between  the  shear  layer 

instability  and  the  feedback  process.  Moreover,  we  note  that  the  low  frequency  peak 
is  located  at  non-dimensional  frequency  Stj^t  =  fD/vg  ss  0.36.  This  is  within  the 
range  0.2  <  Stjet  <  0.5,  where  the  frequency  of  the  preferred  jet  mode  (characteristic 
of  the  largest  scales  on  an  unforced  subsonic  jet)  is  known  to  lie  [12|.  As  we  move 
downstream,  the  amplitudes  for  the  high  frequencies  tend  to  diminish,  as  can  be  ex¬ 
pected.  In  particular.  Figs.  5b,c  indicate  locations  of  the  beginning  and  concluding 
stages,  respectively,  of  a  first  merging  (associated  with  5ti).  Similarly,  Figs.  5d  and  5e 
correspond  to  locations  where  second  (associated  with  Sta)  and  third  (associated  with 
Sta  ss  Stg/S)  mergings  take  place. 


The  detailed  flow  visualization  through  the  sequence  of  panels  in  Fig.  6  allows 
to  determine  the  approximate  pairing  locations,  as  defined  by  the  locations  where  the 
vertical  alignment  of  the  vortices  takes  place.  In  this  way  the  locations  of  the  second 
and  third  pairings  can  be  estimated  from  panels  c)  and  d),  respectively,  at  Z  ss  4.8  cm, 
and  Z  »  9.7  cm.  The  beginning  of  the  first  pairing  can  be  also  observed  in  panel  d), 
at  Z  «  3.3  cm,  where  the  vortices  are  seen  in  the  process  of  rolling  around  each  other. 
There  is  a  close  agreement  between  these  merging  locations  obtained  by  direct  flow 
visualization,  and  the  results  of  the  spectral  analysis  above.  Note,  in  particular,  that 
the  third  merging  actually  occurs  downstream  of  the  end  of  the  potential  core  of  the  jet. 
Because  of  this,  the  effect  of  this  merging  on  the  shear  layer  just  ahead  of  the  trailing 
edge  of  the  nozzle  is  negligible,  as  indicated  by  the  spectral  distribution  in  Fig.  5a. 
This  is  evidence  of  the  role  of  the  jet  diameter,  and  hence  of  the  location  of  the  end  of 
the  potential  core,  in  the  selection  of  the  dominant  low  frequency  mode  of  the  jet. 

Effect  of  Random  Inflow  Perturbations 

We  have  also  simulated  the  effects  of  high  frequency  random  perturbations  in 
the  inflowing  jet  stream.  Such  fluctuations  are  typically  present  in  the  experimental 
conditions  in  the  laboratory,  and  are  due  to  turbulence  and  boundary  layers  in  the 
nozzle.  In  the  calculations  this  is  done  by  replacing  the  streamwise  velocity  U  at 
the  inflow  (i.e.,  for  Z  =  0.0  cm,  R  <  Rg)  with  17(1  4-  p),  where  p  =  p{R,t)  is  the 
perturbating  term.  This  term  is  defined  by  a  sine  Fourier  series  in  the  variable  R,  with 
Rg  the  largest  wavelength: 

1  ^ /? 

P(^.0  =  ~  +  6m),  (6) 

where  ^m(0  is  a  time-dependent  amplitude  defined  by 

dim  i 

«„(0  =  5%=^  +  2,  if  <i<t^  + 

dim  i 


and  by 


otherwise.  In  this  way,  each  added  perturbation  term  varies  between  0.0  and  a  max¬ 
imum  value  of  no  more  than  a  fraction  Pm/m  of  the  inflow  jet  velocity.  The  dura¬ 
tion  6tm,  and  maximum  value  p^  of  the  amplitude,  as  well  as  the  intrinsic  phases 
6m{'^  ^  6m  <  27r)  of  the  terms  in  the  Fourier  series  are  randomly  generated  numbers. 

A  calculation  was  performed  which  restarted  the  program  at  time  t  =  2.52  ms  of 
the  previous  unperturbed  calculation,  but  now  including  inflow  fluctuations.  Here,  the 
series  in  eq.  (6)  was  truncated  to  its  first  four  terms  (M=4),  including  only  fluctuations 
of  wavelength  Rg,  Ro/2,  Rg/^,  and  Rg/A,  with  durations  in  the  range  of  5  —  12.5  pa 
(20  -  50  time  steps).  Due  to  the  changes  in  the  frequency  content  of  the  velocity 
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field,  the  previous  regularity  in  the  flow  pattern  was  lost  in  spite  of  the  relatively 
small  amplitude  of  the  fluctuations.  This  can  be  seen  by  comparison  of  the  spectral 
distributions  obtained  in  this  case  with  those  for  the  unperturbed  case  discussed  above. 
We  compare  Figs.  7a-c,  for  mean  perturbation  levels  of  0.1,  1  and  5  Vi,  respectively,  with 
Fig.  5a.  This  indicates  that  the  peaks  become  broader,  while  new  high  frequency  peaks 
are  now  present.  As  the  fluctuation  level  increases,  the  calculated  spectral  distribution 
begins  to  look  like  the  distributions  obtained  in  the  analysis  of  near  field  pressure 
fluctuations  in  experiments  with  non-excited-jets  [12].  The  experimentally  obtained 
distributions  consist  of  a  broad  peak  tilted  towards  the  low  frequency  side,  and  centered 
on  Stjet  ^  0.37,  in  agreement  with  the  value  of  Stj^t  found  in  the  present  numerical 
simulations. 

Summary  and  Conclusions 

We  have  presented  results  from  numerical  simulations  of  the  evolution  of  the 
Kelvin-Helmholtz  instability  for  a  subsonic,  compressible,  axisymmetric,  spatially- 
evolving  shear  layer.  In  the  absence  of  boundary  layers  and  small-scale  inflow  tur¬ 
bulence,  the  calculations  indicate  that  the  large-scale  development  of  the  unforced  jet 
shear  layer  has  an  underlying  degree  of  organization.  This  is  the  result  of  a  feedback 
mechanism  in  which  the  shear  layer  ahead  of  the  nozzle  edge  is  modulated  by  the  far 
field  induced  by  the  mergings  downstream,  near  the  end  of  the  potential  core  of  the 
jet.  This  is  in  agreement  with  the  experimental  observations  [12]. 

The  dominant  frequencies  in  the  unforced  flow  depend  on  the  two  length-scales  of 
the  jet,  namely,  the  initial  shear  layer  thickness  and  the  jet  diameter.  These  scales  are 
associated  with  the  two  basic  mechanisms  in  the  jet,  the  jet  shear  layer  instability  and 
the  jet  column  instability,  which  determine  the  high  and  low  frequency  natural  modes 
of  the  jet,  respectively.  The  spectral  analysis  of  the  axial  velocity  fluctuations  at  the 
center  of  the  shear  layer,  just  ahead  of  the  nozzle  edge,  shows  a  sequence  of  peaks 
separated  by  a  regular  interval  approximately  equal  to  Stg/A.  These  peaks  correspond 
to  the  subharmonics  of  the  jet  shear  layer  instability  frequency  Sto  and  to  the  non-linear 
interactions  between  them.  The  non-dimensional  frequency  of  the  preferred  jet  mode, 
Stjit  =  (/o/4)i?/vo  ss  0.36,  is  in  good  agreement  with  the  experimental  values  for 
subsonic  jets.  Moreover,  the  value  for  Sto  was  consistent  with  the  predictions  of  linear 
inviscid  instability  theory.  In  addition,  the  distribution  of  the  merging  locations  was 
in  qualitative  agreement  with  that  observed  in  jet-experiments  with  low-level  forcing. 

The  studies  with  random  high  frequency  perturbations  on  the  inflow  velocity  show 
that  such  perturbations  tend  to  break  the  organized  jet  shear  layer,  and  hence  the  tem¬ 
poral  coherence  between  the  structures.  By  generating  more  incoherent  mergings  the 
fluctuations  tend  to  destroy  the  regularity  of  the  idealized,  unforced  two-dimensional 
case. 
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Fig.  2:  Flow  configuration  for  the  axisymmelric  jet 
simulation.  The  ’x'  indicates  the  location  of  the  trailing  edge 
of  the  nozzle. 
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Fig.  5:  Spectral  amplitude  at  R  =  D/2,  and  van 
ous  axial  locatioru,  as  a  function  of  Strouhal  frequency  St 
a)  Z  =  2.39  cm  (z  =  0.40);  h)  Z  =  2.94  cm  (z  =  0.71) 
c)  Z  =  3.54  cm  (z  =  1.04);  d)  Z  =  5.24  cm  (z  =  1.98) 
«)  Z  =  9.44  cm  (z  =  4.32). 
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